rm(list=ls(all=TRUE))
library(arm)#for using "attach.all", "lmer", and "sim" functions
library(foreign)#for reading in Stata datasets

#create function for standard error of a mean
mean.se <- function (x, na.rm) {stderr <- sd(x, na.rm = TRUE)
    return(stderr/sqrt(sum(!is.na(x))))}

###################
#get data on proportion of REPs in each circuit-year
partisan.comp <- read.csv("partisan.composition.csv")
year <-partisan.comp$year
circuit <- partisan.comp$circuit
circuit.size <- partisan.comp$active.judges.circyear
dem.proportion <- partisan.comp$dem.prop.active.circyear #proportion of active judges
gop.proportion <- 1 -dem.proportion
dem.circuit <- ifelse(partisan.comp$dem.prop.active.circyear > .5,1,0)
gop.circuit <-  ifelse(partisan.comp$dem.prop.active.circyear < .5,1,0)
split.circuit <-  ifelse(partisan.comp$dem.prop.active.circyear ==.5,1,0)
partisan.comp.merge <- data.frame(year, circuit, circuit.size, #create small dataset for merging below
	dem.proportion, gop.proportion, dem.circuit, gop.circuit, split.circuit) 


#JUDGE-LEVEL DATA

judge.data <- read.dta("data_for_analysis_judge_level.dta", convert.underscore=T)
attach.all(judge.data, overwrite = TRUE)
#keep only needed variables
judge.data.keep <- data.frame(year, circuit, caseid, casename, judge.indicator, citation, issue, lowervote,
  panelvote, dissent, direct.from.agency, issue.name, j1name, j1vote, j1code, j1jcs, j1dem, 
  j2name, j2vote, j2code, j2jcs, j2dem,  j3name, j3vote, j3code, j3jcs, j3dem)
dim(judge.data.keep)
#merge partisan composition data
judge.data.keep <- merge(judge.data.keep, partisan.comp.merge)
dim(judge.data.keep)

#code panels based on j2 and j3
  #use "coll" as short for "colleagues"
judge.data.keep$dd.coll <- ifelse(judge.data.keep$j2dem + judge.data.keep$j3dem == 2,1,0)
judge.data.keep$dr.coll <- ifelse(judge.data.keep$j2dem + judge.data.keep$j3dem == 1,1,0)
judge.data.keep$rr.coll <- ifelse(judge.data.keep$j2dem + judge.data.keep$j3dem == 0,1,0)

#define panel majorties in terms of J1-- that is only code cases where J1 is in majority
judge.data.keep$dem.majority <- ifelse((judge.data.keep$j2dem ==1 | judge.data.keep$j3dem == 1) & judge.data.keep$j1dem == 1, 1,0)    #define majority democratic panels 
judge.data.keep$gop.majority <- ifelse((judge.data.keep$j2dem ==0 | judge.data.keep$j3dem == 0)& judge.data.keep$j1dem == 0,1,0)#define GOP-majority panel
#define majority panels, in terms of j1
#for both dems and gop separately, use nearest neighbor matching to check whethe we need to create balanced data set,
  #on majority panels
  #where counterj is addition of a minority member to panel
  #that is, for DEMS, adding a Republican instead of a third Democrat
  #that is, for GOP, adding a Democrat instead of a third Republican
judge.data.keep$gop.counterj <-  ifelse((judge.data.keep$j2dem ==0 | judge.data.keep$j3dem == 0) & judge.data.keep$dem.majority ==1, 1,0)
judge.data.keep$dem.counterj <-  ifelse((judge.data.keep$j2dem ==1 | judge.data.keep$j3dem == 1) & judge.data.keep$gop.majority == 1, 1,0)
detach()
attach.all(judge.data.keep, overwrite = "T")

# issue.name <- judge.data.keep$issue.name
# unique.issue <- sort(unique(issue.name))
# case.name <- judge.data.keep$casename
# unique.case <- sort(unique(case.name))
# J <- length(unique.issue)

#get differences in voting rates between dems and rep appointeees
dem.rate <- mean(j1vote[j1dem==1])
gop.rate <- mean(j1vote[j1dem==0])
options(digits=2)
print(c(dem.rate, gop.rate))

#pooling all judges, get votes based on different colleagues

all.dd.rate <- mean(j1vote[dd.coll==1])
all.dd.se <-  mean.se(j1vote[dd.coll==1])
all.dr.rate <- mean(j1vote[dr.coll==1])
all.dr.se <-  mean.se(j1vote[dr.coll==1])
all.rr.rate <- mean(j1vote[rr.coll==1])
all.rr.se <-  mean.se(j1vote[rr.coll==1])

dem.dd.rate <- mean(j1vote[dd.coll==1 & j1dem == 1])
dem.dd.se <-  mean.se(j1vote[dd.coll==1 & j1dem == 1])
dem.dr.rate <- mean(j1vote[dr.coll==1 & j1dem == 1])
dem.dr.se <-  mean.se(j1vote[dr.coll==1 & j1dem == 1])
dem.rr.rate <- mean(j1vote[rr.coll==1 & j1dem == 1])
dem.rr.se <-  mean.se(j1vote[rr.coll==1 & j1dem == 1])

gop.dd.rate <- mean(j1vote[dd.coll==1 & j1dem == 0])
gop.dd.se <-  mean.se(j1vote[dd.coll==1 & j1dem == 0])
gop.dr.rate <- mean(j1vote[dr.coll==1 & j1dem == 0])
gop.dr.se <-  mean.se(j1vote[dr.coll==1 & j1dem == 0])
gop.rr.rate <- mean(j1vote[rr.coll==1 & j1dem == 0])
gop.rr.se <-  mean.se(j1vote[rr.coll==1 & j1dem == 0])

#RUN SCRIPT TO CREATE FIGURE 1
source("figure_1_script.R")

#get percent of affirmances
mean(j1vote == lowervote)#64%


#RUN SCRIPT TO MAKE FIGURE 3a
source("figure_3a_script.R")

#RUN SCRIPT TO MAKE FIGURE 3b
source("figure_3b_script.R")


################################################

#MULTILEVEL MODELS USED IN TABLE 1
##for Dems and GOP separately, estimate the counterj effect of switching from a unified panel 
  #to panel with 1 other member of the other party
  #define "counterj" as addition of judge from other party
    #so gop.counterj == addition of a gop judge compared to Democrat

sims <- 1000#number of replications in simulations

##################################

#A) 1ST COLUMN:  FOR GOP MAJORITIES
gop.majority.data <- judge.data.keep[judge.data.keep$gop.majority == 1,]
detach()
attach.all(gop.majority.data, overwrite = T)

z.j1name <- as.factor(as.character(j1name)) #rename and code judges to make it easier to pull out random effects
mlm.gop <- glmer(j1vote ~ j1jcs  + lowervote  + dem.counterj +  
  dem.circuit + dem.circuit:dem.counterj  +  (1|issue.name) + (1|circuit) + (1|year) + (1|z.j1name),
   family=binomial(link="logit"))
display(mlm.gop)
ranef(mlm.gop)
se.ranef(mlm.gop)
pred.mlm.gop <- ifelse(fitted(mlm.gop) >=.5, 1,0)#get prediction from model
modal.cat <- ifelse(mean(j1vote) > .5, mean(j1vote), 1- mean(j1vote))#percent in modal category; flips if mean < .5
correct.mlm.gop <- pred.mlm.gop == j1vote #obsevervations prediced correctly
correct.percent.mlm.gop <- sum(correct.mlm.gop)/length(j1vote)#percent correct
prop.reduct.mlm.gop <- 100*((100*correct.percent.mlm.gop - 100*modal.cat)/(100-100*modal.cat)) #PRE
correct.percent.mlm.gop
prop.reduct.mlm.gop
sim.mlm.gop <-  sim(mlm.gop, n.sims=sims)#simulate model 1000 times

#get total effect of dem.counterj in dem circuits
betas.dem.counterj <- sim.mlm.gop$fixef[,"dem.counterj"]
betas.dem.interaction <- sim.mlm.gop$fixef[,"dem.counterj:dem.circuit"]
total.dem.counterj <- quantile(betas.dem.counterj + betas.dem.interaction, c(.025,.5,.975))
total.dem.counterj
#get total effect of dem circuit when dem.counterj ==1
betas.dem.circuit <- sim.mlm.gop$fixef[,"dem.circuit"]
total.dem.circuit <- quantile(betas.dem.circuit + betas.dem.interaction, c(.025,.5,.975))
total.dem.circuit

##################################

#B) 2nd COLUMN:  FOR DEM MAJORITIES
#create separate data frame to make it easier to pull out random effects below
dem.majority.data <- judge.data.keep[judge.data.keep$dem.majority == 1,]
detach()
attach.all(dem.majority.data, overwrite = T)
dim(dem.majority.data)

z.j1name <- as.factor(as.character(j1name)) #rename and code judges to make it easier to pull out random effects
mlm.dems <- glmer(j1vote ~ j1jcs  + lowervote + gop.circuit + gop.counterj +
    gop.circuit:gop.counterj
    + (1|issue.name) + (1|circuit) + (1|year) + (1|z.j1name)
   , family=binomial(link="logit"))
display(mlm.dems)
ranef(mlm.dems)
se.ranef(mlm.dems)
pred.mlm.dems <- ifelse(fitted(mlm.dems) >=.5, 1,0)#get prediction from model
modal.cat <- ifelse(mean(j1vote) > .5, mean(j1vote), 1- mean(j1vote))#percent in modal category; flips if mean < .5
correct.mlm.dems <- pred.mlm.dems == j1vote#obsevervations prediced correctly
correct.percent.mlm.dems <- sum(correct.mlm.dems)/length(j1vote)#percent correct
prop.reduct.mlm.dems <- 100*((100*correct.percent.mlm.dems - 100*modal.cat)/(100-100*modal.cat))#PRE
aic.mlm.dems <- AIC(logLik(mlm.dems)) 
modal.cat
correct.percent.mlm.dems
prop.reduct.mlm.dems
sim.mlm.dems <-  sim(mlm.dems, n.sims=sims)

#get total effect of gop.counterj in gop circuits
betas.gop.counterj <- sim.mlm.dems$fixef[,"gop.counterj"]
betas.gop.interaction <- sim.mlm.dems$fixef[,"gop.circuit:gop.counterj"]
total.gop.counterj <- quantile(betas.gop.counterj + betas.gop.interaction, c(.025,.5,.975))
total.gop.counterj
 #get total effect of gop circuit when gop.counterj ==1
betas.gop.circuit <- sim.mlm.dems$fixef[,"gop.circuit"]
total.gop.circuit <- quantile(betas.gop.circuit + betas.gop.interaction, c(.025,.5,.975))
total.gop.circuit 

####
#CREATE ESTIMATES FOR FIGURE 3:
#A) Republican majorities: COMPUTE PREDICTED PROBS AND AVERAGE PREDICTIVE DIFFERENCE OF DEM counterj on GOP JUDGES 
	#IN DEM CIRCUITS AND GOP CIRCUITS
yes <- 1 #indicator for counterj
no <- 0
gop.judges.apd.dem.circuit.vector <- rep(NA, sims)
gop.judges.apd.gop.circuit.vector <- rep(NA, sims)
gop.judges.prob.dem.circuit.mixed.vector <- rep(NA, sims)
gop.judges.prob.dem.circuit.unified.vector <- rep(NA, sims)
gop.judges.prob.gop.circuit.mixed.vector <- rep(NA, sims)
gop.judges.prob.gop.circuit.unified.vector <- rep(NA, sims)

  # note that sim returns:
#    $fixef (fixed effects)
#    $j1name
#    $year
#    $circuit
#    $issue.name
    #define ranef for sims


for (i in 1:sims){#loop over simulations
    b <- sim.mlm.gop$fixef[i,]#grab "fixed effects"
        year.re <- sim.mlm.gop$year[i,,][as.factor(year)]
        judge.re <- sim.mlm.gop$z.j1name[i,,][as.factor(z.j1name)]
        circuit.re <- sim.mlm.gop$circuit[i,,][as.factor(circuit)]
        issue.re <- sim.mlm.gop$issue[i,,][as.factor(issue)]      
    #get prediction when dem.counterj == 1 and gop.circuit = 1
    gop.judges.prob.gop.circuit.mixed <-  invlogit(b["(Intercept)"] + b["j1jcs"]*j1jcs + 
    b["lowervote"]*lowervote +   b["dem.counterj"]*yes   
    + year.re + issue.re + circuit.re +  judge.re) #   add random effects
    #get prediction when dem.counterj == 0 and gop.circuit = 1
    gop.judges.prob.gop.circuit.unified <-  invlogit(b["(Intercept)"] + b["j1jcs"]*j1jcs + b["lowervote"]*lowervote +
    b["dem.counterj"]*no    + year.re + issue.re + circuit.re + judge.re) 
    delta.gop.circuit <- gop.judges.prob.gop.circuit.mixed - gop.judges.prob.gop.circuit.unified 
    #################################
    #NOW DO IN DEM CIRCUITS (add in interaction)
    #dem.counterj == 1
    gop.judges.prob.dem.circuit.mixed <-  invlogit(b["(Intercept)"] + b["j1jcs"]*j1jcs + b["lowervote"]*lowervote +
    b["dem.counterj"]*yes   + b["dem.circuit"]*yes + b["dem.counterj:dem.circuit"]*yes
    + year.re + issue.re + circuit.re + judge.re) 
    #dem.counterj == 0 
    gop.judges.prob.dem.circuit.unified <-   invlogit(b["(Intercept)"] + b["j1jcs"]*j1jcs + 
    b["lowervote"]*lowervote +  b["dem.counterj"]*no   + b["dem.circuit"]*yes 
    + b["dem.counterj:dem.circuit"]*no + year.re + issue.re + circuit.re +  judge.re) 
    delta.dem.circuit <- gop.judges.prob.dem.circuit.mixed - gop.judges.prob.dem.circuit.unified
    gop.judges.prob.dem.circuit.unified.vector[i] <- mean(gop.judges.prob.dem.circuit.unified)
    gop.judges.prob.dem.circuit.mixed.vector[i] <- mean(gop.judges.prob.dem.circuit.mixed)
    gop.judges.apd.dem.circuit.vector[i] <- mean(delta.dem.circuit)
    gop.judges.prob.gop.circuit.unified.vector[i] <- mean(gop.judges.prob.gop.circuit.unified)
    gop.judges.prob.gop.circuit.mixed.vector[i] <- mean(gop.judges.prob.gop.circuit.mixed)
    gop.judges.apd.gop.circuit.vector[i] <- mean(delta.gop.circuit)
}

#average differences
quantile(gop.judges.apd.dem.circuit.vector, c(.025, .5, .975))
quantile(gop.judges.apd.gop.circuit.vector, c(.025, .5, .975))
#predicted probabilities
quantile(gop.judges.prob.dem.circuit.unified.vector, c(.025, .5, .975))
quantile(gop.judges.prob.dem.circuit.mixed.vector, c(.025, .5, .975))
quantile(gop.judges.prob.gop.circuit.unified.vector, c(.025, .5, .975))
quantile(gop.judges.prob.gop.circuit.mixed.vector, c(.025, .5, .975))

##################

#B) DEMOCRATIC Majorities: COMPUTE PREDICTED PROBS AND AVERAGE PREDICTIVE DIFFERENCE 
	#OF GOP counterj on Dem majorities IN DEM CIRCUITS AND GOP CIRCUITS

yes <- 1 #indicator for counterj
no <- 0
#create empty vectors to store results
dem.judges.apd.dem.circuit.vector <- rep(NA, sims)
dem.judges.apd.gop.circuit.vector <- rep(NA, sims)
dem.judges.prob.dem.circuit.mixed.vector <- rep(NA, sims)
dem.judges.prob.dem.circuit.unified.vector <- rep(NA, sims)
dem.judges.prob.gop.circuit.mixed.vector <- rep(NA, sims)
dem.judges.prob.gop.circuit.unified.vector <- rep(NA, sims)

for (i in 1:sims){#loop over simulations
    b <- sim.mlm.dems$fixef[i,]#grab "fixed effects"
    #need to pull out proper ranefs, by each run of simulation, then merge to data-length
    year.re <- sim.mlm.dems$year[i,,][as.factor(year)]
    judge.re <- sim.mlm.dems$z.j1name[i,,][as.factor(z.j1name)]
    circuit.re <- sim.mlm.dems$circuit[i,,][as.factor(circuit)]
    issue.re <- sim.mlm.dems$issue[i,,][as.factor(issue)]

    #get prediction when gop.counterj == 1 and dem.circuit = 1
    dem.judges.prob.dem.circuit.mixed <-  invlogit(b["(Intercept)"] + b["j1jcs"]*j1jcs +
    	b["lowervote"]*lowervote +  b["gop.counterj"]*yes   
    	+ year.re + issue.re + circuit.re +  judge.re) #   add random effects
    #get prediction when gop.counterj == 0 and dem.circuit = 1
    dem.judges.prob.dem.circuit.unified <-  invlogit(b["(Intercept)"] + b["j1jcs"]*j1jcs + b["lowervote"]*lowervote +
    	b["gop.counterj"]*no  + year.re + issue.re + circuit.re + judge.re) 
    	delta.dem.circuit <- dem.judges.prob.dem.circuit.mixed - dem.judges.prob.dem.circuit.unified 
    #################################
    #NOW DO IN GOP CIRCUITS (add in interaction)
    #gop.counterj == 1
    dem.judges.prob.gop.circuit.mixed <-  invlogit(b["(Intercept)"] + b["j1jcs"]*j1jcs + b["lowervote"]*lowervote +
    b["gop.counterj"]*yes   + b["gop.circuit"]*yes + b["gop.circuit:gop.counterj"]*yes
    + year.re + issue.re + circuit.re +  judge.re) 
    #gop.counterj == 0 
    dem.judges.prob.gop.circuit.unified <-   invlogit(b["(Intercept)"] + b["j1jcs"]*j1jcs + 
    b["lowervote"]*lowervote +  b["gop.counterj"]*no   + b["gop.circuit"]*yes 
    + b["gop.circuit:gop.counterj"]*no + year.re + issue.re + circuit.re +  judge.re) 
    delta.gop.circuit <- dem.judges.prob.gop.circuit.mixed - dem.judges.prob.gop.circuit.unified
    
    dem.judges.prob.dem.circuit.unified.vector[i] <- mean(dem.judges.prob.dem.circuit.unified)
    dem.judges.prob.dem.circuit.mixed.vector[i] <- mean(dem.judges.prob.dem.circuit.mixed)
    dem.judges.apd.dem.circuit.vector[i] <- mean(delta.dem.circuit)
    dem.judges.prob.gop.circuit.unified.vector[i] <- mean(dem.judges.prob.gop.circuit.unified)
    dem.judges.prob.gop.circuit.mixed.vector[i] <- mean(dem.judges.prob.gop.circuit.mixed)
    dem.judges.apd.gop.circuit.vector[i] <- mean(delta.gop.circuit)
    }

#average differences
quantile(dem.judges.apd.dem.circuit.vector, c(.025, .5, .975))
quantile(dem.judges.apd.gop.circuit.vector, c(.025, .5, .975))
#predicted probabilities
quantile(dem.judges.prob.dem.circuit.unified.vector, c(.025, .5, .975))
quantile(dem.judges.prob.dem.circuit.mixed.vector, c(.025, .5, .975))
quantile(dem.judges.prob.gop.circuit.unified.vector, c(.025, .5, .975))
quantile(dem.judges.prob.gop.circuit.mixed.vector, c(.025, .5, .975))

##########
#RUN SCRIPT TO CREATE FIGURE 4a
source("figure_4a_script.R")

####################################
####################################
#CASE-LEVEL DATA 
###############################################
case.data <- read.dta("data_for_analysis_case_level.dta", convert.underscore=T)
#merge partisan composition data
dim(case.data)
case.data <- merge(case.data, partisan.comp.merge)
dim(case.data)
detach()
attach.all(case.data, overwrite = TRUE)

#########
#CODE, AGGREGATE PANEL RATES, FIRST FOR ALL CASES, THEN BY ISSUE
rrr.panel <- ifelse(j1dem == 0 & j2dem == 0 & j3dem == 0,1,0)
rrd.panel <- ifelse(j1dem+j2dem+j3dem == 1,1,0)
ddr.panel <- ifelse(j1dem+j2dem+j3dem == 2,1,0)
ddd.panel <- ifelse(j1dem == 1 & j2dem == 1 & j3dem == 1,1,0)
all.panel <- ifelse(rrr.panel == 1, "rrr", 
    ifelse(rrd.panel==1, "rrd",
        ifelse(ddr.panel == 1, "ddr",
            ifelse(ddd.panel == 1, "ddd", NA))))
table(rrr.panel)
table(rrd.panel)
table(ddr.panel)
table(ddd.panel)

rrr.rate <- mean(panelvote[rrr.panel==1])
rrr.se <- mean.se(panelvote[rrr.panel==1])
rrd.rate <- mean(panelvote[rrd.panel==1])
rrd.se <- mean.se(panelvote[rrd.panel==1])
ddr.rate <- mean(panelvote[ddr.panel==1])
ddr.se <- mean.se(panelvote[ddr.panel==1])
ddd.rate <- mean(panelvote[ddd.panel==1])
ddd.se <- mean.se(panelvote[ddd.panel==1])

#break down by circuits 

rrr.dem.circ.rate <- mean(panelvote[rrr.panel==1 & dem.circuit == 1])
rrr.dem.circ.se <- mean.se(panelvote[rrr.panel==1 & dem.circuit == 1])
rrd.dem.circ.rate <- mean(panelvote[rrd.panel==1 & dem.circuit == 1])
rrd.dem.circ.se <- mean.se(panelvote[rrd.panel==1 & dem.circuit == 1])
ddr.dem.circ.rate <- mean(panelvote[ddr.panel==1 & dem.circuit == 1])
ddr.dem.circ.se <- mean.se(panelvote[ddr.panel==1 & dem.circuit == 1])
ddd.dem.circ.rate <- mean(panelvote[ddd.panel==1 & dem.circuit == 1])
ddd.dem.circ.se <- mean.se(panelvote[ddd.panel==1 & dem.circuit == 1])

rrr.gop.circ.rate <- mean(panelvote[rrr.panel==1 & gop.circuit == 1])
rrr.gop.circ.se <- mean.se(panelvote[rrr.panel==1 & gop.circuit == 1])
rrd.gop.circ.rate <- mean(panelvote[rrd.panel==1 & gop.circuit == 1])
rrd.gop.circ.se <- mean.se(panelvote[rrd.panel==1 & gop.circuit == 1])
ddr.gop.circ.rate <- mean(panelvote[ddr.panel==1 & gop.circuit == 1])
ddr.gop.circ.se <- mean.se(panelvote[ddr.panel==1 & gop.circuit == 1])
ddd.gop.circ.rate <- mean(panelvote[ddd.panel==1 & gop.circuit == 1])
ddd.gop.circ.se <- mean.se(panelvote[ddd.panel==1 & gop.circuit == 1])

mean.vector <- c(ddd.rate, ddr.rate,rrd.rate,  rrr.rate)
se.vector <-  c(ddd.se, ddr.se, rrd.se, rrr.se)

mean.dem.circ.vector <- c(ddd.dem.circ.rate, ddr.dem.circ.rate,rrd.dem.circ.rate,  rrr.dem.circ.rate)
se.dem.circ.vector <-  c(ddd.dem.circ.se, ddr.dem.circ.se,rrd.dem.circ.se,  rrr.dem.circ.se)
mean.gop.circ.vector <- c(ddd.gop.circ.rate, ddr.gop.circ.rate,rrd.gop.circ.rate,  rrr.gop.circ.rate)
se.gop.circ.vector <-  c(ddd.gop.circ.se, ddr.gop.circ.se,rrd.gop.circ.se,  rrr.gop.circ.se)
#add hypothetical for pure majoritarian control world
mean.vector.hyp <- c(mean.vector[1], mean.vector[1], mean.vector[4],mean.vector[4])
mean.dem.circ.vector.hyp <- c(mean.dem.circ.vector[1],mean.dem.circ.vector[1],  mean.dem.circ.vector[4],mean.dem.circ.vector[4])
mean.gop.circ.vector.hyp <- c(mean.gop.circ.vector[1],mean.gop.circ.vector[1],  mean.gop.circ.vector[4],mean.gop.circ.vector[4])

#Make Figure 4

source("figure_4b_script.R")


##############


#####################################################################
#ROBUSTNESS CHECKS (PRESENTED IN TABLE A-2)

##########
#FOR GOP MAJORITIES
detach()
attach.all(gop.majority.data, overwrite = T)
z.j1name <- as.factor(as.character(j1name)) #rename and code judges to make it easier to pull out random effects

#original model
mlm.gop <- glmer(j1vote ~ j1jcs  + lowervote  +   dem.circuit 
		+ dem.counterj + dem.circuit:dem.counterj  +  (1|issue.name) + (1|circuit) + (1|year) + (1|z.j1name),
   family=binomial(link="logit"))
display(mlm.gop, detail = T)
pred.mlm.gop <- ifelse(fitted(mlm.gop) >=.5, 1,0)#get prediction from model
modal.cat <- ifelse(mean(j1vote) > .5, mean(j1vote), 1- mean(j1vote))#flips if mean < .5
correct.mlm.gop <- pred.mlm.gop == j1vote
correct.percent.mlm.gop <- sum(correct.mlm.gop)/length(j1vote)
prop.reduct.mlm.gop <- 100*((100*correct.percent.mlm.gop - 100*modal.cat)/(100-100*modal.cat))
correct.percent.mlm.gop
prop.reduct.mlm.gop
#aic.mlm.gop <- AIC(logLik(mlm.gop)) 
sim.mlm.gop <-  sim(mlm.gop, n.sims=sims)
betas.dem.counterj <- sim.mlm.gop$fixef[,"dem.counterj"]
betas.dem.interaction <- sim.mlm.gop$fixef[,"dem.circuit:dem.counterj"]
total.dem.counterj <- quantile(betas.dem.counterj + betas.dem.interaction, c(.025,.5,.975))
total.dem.counterj

#model 2: dropping Re's for year and judge
mlm.gop.2 <- glmer(j1vote ~ j1jcs  + lowervote  +  dem.circuit + dem.counterj +  
  dem.circuit:dem.counterj  +  (1|issue.name) + (1|circuit),
   family=binomial(link="logit"))
display(mlm.gop.2, detail = T)
pred.mlm.gop.2 <- ifelse(fitted(mlm.gop.2) >=.5, 1,0)
modal.cat <- ifelse(mean(j1vote) > .5, mean(j1vote), 1- mean(j1vote))#flips if mean < .5
correct.mlm.gop.2 <- pred.mlm.gop.2 == j1vote
correct.percent.mlm.gop.2 <- sum(correct.mlm.gop.2)/length(j1vote)
prop.reduct.mlm.gop.2 <- 100*((100*correct.percent.mlm.gop.2 - 100*modal.cat)/(100-100*modal.cat))
correct.percent.mlm.gop.2
prop.reduct.mlm.gop.2

#regular logit, with fixed effects for issues and years

logit.gop.1 <- glm(j1vote ~ j1jcs  + lowervote + dem.circuit + dem.counterj +
	dem.circuit:dem.counterj + as.factor(issue) + as.factor(year), family=binomial(link="logit"))
display(logit.gop.1, detail = T)
pred.logit.gop.1 <- ifelse(fitted(logit.gop.1) >=.5, 1,0)
modal.cat <- ifelse(mean(j1vote) > .5, mean(j1vote), 1- mean(j1vote))#flips if mean < .5
correct.logit.gop.1 <- pred.logit.gop.1 == j1vote
correct.percent.logit.gop.1 <- sum(correct.logit.gop.1)/length(j1vote)
prop.reduct.logit.gop.1 <- 100*((100*correct.percent.logit.gop.1 - 100*modal.cat)/(100-100*modal.cat))
correct.percent.logit.gop.1
prop.reduct.logit.gop.1


#use GEE, with cases as clusters, with fixed effects for issues  and years
gop.majority.data.sorted <- gop.majority.data[order(gop.majority.data$caseid),] #sort by caseid
detach()
attach.all(gop.majority.data.sorted)
#with years and issues
library(geepack)
gee.gop.1 <- geeglm(j1vote ~ j1jcs  + lowervote + dem.circuit + dem.counterj +
	dem.circuit:dem.counterj + as.factor(issue) + as.factor(year), data=gop.majority.data.sorted,
	id =gop.majority.data.sorted$caseid, corstr = "exchangeable", family=binomial)
summary(gee.gop.1)
pred.gee.gop.1 <- ifelse(fitted(gee.gop.1) >=.5, 1,0)
modal.cat <- ifelse(mean(j1vote) > .5, mean(j1vote), 1 - mean(j1vote))#flips if mean < .5
correct.gee.gop.1 <- pred.gee.gop.1 == j1vote
correct.percent.gee.gop.1 <- sum(correct.gee.gop.1)/length(j1vote)
prop.reduct.gee.gop.1 <- 100*((100*correct.percent.gee.gop.1 - 100*modal.cat)/(100-100*modal.cat))
modal.cat
correct.percent.gee.gop.1
prop.reduct.gee.gop.1

############################################################
############################################################
#FOR DEM MAJORITIES
detach()
attach.all(dem.majority.data, overwrite = T)
dim(dem.majority.data)

#original model
z.j1name <- as.factor(as.character(j1name)) #rename and code judges to make it easier to pull out random effects
mlm.dems <- glmer(j1vote ~ j1jcs  + lowervote + gop.circuit + gop.counterj +
    gop.circuit:gop.counterj
    + (1|issue.name) + (1|circuit) + (1|year) + (1|z.j1name)
   , family=binomial(link="logit"))
display(mlm.dems, detail = T)
pred.mlm.dems <- ifelse(fitted(mlm.dems) >=.5, 1,0)
modal.cat <- ifelse(mean(j1vote) > .5, mean(j1vote), 1- mean(j1vote))#flips if mean < .5
correct.mlm.dems <- pred.mlm.dems == j1vote
correct.percent.mlm.dems <- sum(correct.mlm.dems)/length(j1vote)
prop.reduct.mlm.dems <- 100*((100*correct.percent.mlm.dems - 100*modal.cat)/(100-100*modal.cat))
aic.mlm.dems <- AIC(logLik(mlm.dems)) 
aic.mlm.dems
modal.cat
correct.percent.mlm.dems
prop.reduct.mlm.dems

#dropping Re's for year and judge
mlm.dems.2 <- glmer(j1vote ~ j1jcs  + lowervote + gop.circuit + gop.counterj +
    gop.circuit:gop.counterj
    + (1|issue.name) + (1|circuit) 
   , family=binomial(link="logit"))
display(mlm.dems.2, detail=T)
pred.mlm.dems.2 <- ifelse(fitted(mlm.dems.2) >=.5, 1,0)
modal.cat <- ifelse(mean(j1vote) > .5, mean(j1vote), 1- mean(j1vote))#flips if mean < .5
correct.mlm.dems.2 <- pred.mlm.dems.2 == j1vote
correct.percent.mlm.dems.2 <- sum(correct.mlm.dems.2)/length(j1vote)
prop.reduct.mlm.dems.2 <- 100*((100*correct.percent.mlm.dems.2 - 100*modal.cat)/(100-100*modal.cat))
aic.mlm.dems.2 <- AIC(logLik(mlm.dems.2)) 
aic.mlm.dems.2
modal.cat
correct.percent.mlm.dems.2
prop.reduct.mlm.dems.2

#regular logit, with fixed issue effects
logit.dem.1 <- glm(j1vote ~ j1jcs  + lowervote + gop.circuit + gop.counterj +
    gop.circuit:gop.counterj + as.factor(issue) + as.factor(year), family=binomial(link="logit"))
display(logit.dem.1, detail = T)
pred.logit.dem.1 <- ifelse(fitted(logit.dem.1) >=.5, 1,0)
modal.cat <- ifelse(mean(j1vote) > .5, mean(j1vote), 1- mean(j1vote))#flips if mean < .5
correct.logit.dem.1 <- pred.logit.dem.1 == j1vote
correct.percent.logit.dem.1 <- sum(correct.logit.dem.1)/length(j1vote)
prop.reduct.logit.dem.1 <- 100*((100*correct.percent.logit.dem.1 - 100*modal.cat)/(100-100*modal.cat))
correct.percent.logit.dem.1
prop.reduct.logit.dem.1

#use GEE, with cases as clusters, with fixed effects for issues 
dem.majority.data.sorted <- dem.majority.data[order(dem.majority.data$caseid),]
detach()
attach.all(dem.majority.data.sorted)
library(geepack)
gee.dems.1 <- geeglm(j1vote ~ j1jcs  + lowervote + gop.circuit + gop.counterj +
	gop.circuit:gop.counterj + as.factor(issue) + as.factor(year), data=dem.majority.data.sorted,
	id =dem.majority.data.sorted$caseid, corstr = "exchangeable", family=binomial)
summary(gee.dems.1)
pred.gee.dems.1 <- ifelse(fitted(gee.dems.1) >=.5, 1,0)
modal.cat <- ifelse(mean(j1vote) > .5, mean(j1vote), 1- mean(j1vote))#flips if mean < .5
correct.gee.dems.1 <- pred.gee.dems.1 == j1vote
correct.percent.gee.dems.1 <- sum(correct.gee.dems.1)/length(j1vote)
prop.reduct.gee.dems.1 <- 100*((100*correct.percent.gee.dems.1 - 100*modal.cat)/(100-100*modal.cat))
correct.percent.gee.dems.1
prop.reduct.gee.dems.1

